{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division, absolute_import\n",
    "  \n",
    "import GPy\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import safeopt\n",
    "\n",
    "mpl.rcParams['figure.figsize'] = (20.0, 10.0)\n",
    "mpl.rcParams['font.size'] = 20\n",
    "mpl.rcParams['lines.markersize'] = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define a kernel and function\n",
    "\n",
    "Here we define a kernel. The function is drawn at random from the GP and is corrupted my Gaussian noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Measurement noise\n",
    "noise_var = 0.05 ** 2\n",
    "\n",
    "# Bounds on the inputs variable\n",
    "bounds = [(-5., 5.), (-5., 5.)]\n",
    "parameter_set = safeopt.linearly_spaced_combinations([bounds[0]], 1000)\n",
    "\n",
    "# Define Kernel\n",
    "# works on the first column of X, index=0\n",
    "k_parameters = GPy.kern.RBF(input_dim=1, variance=2., lengthscale=1.0, active_dims=[0])\n",
    "# works on the second column of X, index=1\n",
    "k_context = GPy.kern.RBF(input_dim=1, variance=2., lengthscale=1.0, active_dims=[1], name='context')\n",
    "kernel = k_parameters * k_context\n",
    "\n",
    "# set of parameters\n",
    "num_contexts = 1\n",
    "\n",
    "# Initial safe point\n",
    "x0 = np.array([[0]])\n",
    "\n",
    "# Generate function with safe initial point at x=0\n",
    "def sample_safe_fun(context=0):\n",
    "    context = np.atleast_2d(context)\n",
    "    while True:\n",
    "        # Joint function over parameters and contexts\n",
    "        sampled_fun = safeopt.sample_gp_function(kernel.copy(), bounds, noise_var, 10)\n",
    "        \n",
    "        if sampled_fun(np.hstack([x0, context]), noise=False) > 0.5:\n",
    "            break\n",
    "        \n",
    "    return sampled_fun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot(context):\n",
    "    # Make points transparent when they belong to a different context\n",
    "    context = np.atleast_2d(context)\n",
    "    opt.context = context\n",
    "    \n",
    "    relevance = opt.gp.kern.context.K(np.hstack([[[0]], context]), opt.gp.X)\n",
    "    relevance /= opt.gp.kern.context.variance\n",
    "    relevance = np.exp(100 * relevance) / np.exp(100)\n",
    "    relevance[relevance < 0.25] = 0.25\n",
    "    point_color = np.zeros((opt.gp.X.shape[0], 4))\n",
    "    point_color[:, 3] = relevance\n",
    "            \n",
    "    # Plot GP\n",
    "    opt.plot(n_samples=1000, point_color=point_color)\n",
    "    \n",
    "    # Plot the true function\n",
    "    data = np.concatenate((parameter_set, np.broadcast_to(context, (parameter_set.shape[0], context.shape[1]))), axis=1)\n",
    "    plt.plot(parameter_set, fun(data, noise=False), color='C2', alpha=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "context = np.array([[0]])\n",
    "\n",
    "# Define the objective function\n",
    "fun = sample_safe_fun(context)\n",
    "\n",
    "# The statistical model of our objective function\n",
    "x = np.hstack([x0, context])\n",
    "gp = GPy.models.GPRegression(x, fun(x), kernel, noise_var=noise_var)\n",
    "\n",
    "# The optimization routine\n",
    "opt = safeopt.SafeOpt(gp, parameter_set, 0., num_contexts=1, threshold=0.5)\n",
    "opt.context = context\n",
    "\n",
    "plot(context)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain next query point\n",
    "x_next = opt.optimize(context)\n",
    "\n",
    "# Get a measurement from the real system\n",
    "y_meas = fun(np.hstack((x_next[None], context)))\n",
    "\n",
    "# Add this to the GP model\n",
    "opt.add_new_data_point(x_next, y_meas, context=context)\n",
    "\n",
    "plot(context=context)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "context = np.array([[0.1]])\n",
    "plot(context)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain next query point\n",
    "x_next = opt.optimize(context)\n",
    "\n",
    "# Get a measurement from the real system\n",
    "y_meas = fun(np.hstack((x_next[None], context)))\n",
    "\n",
    "# Add this to the GP model\n",
    "opt.add_new_data_point(x_next, y_meas, context=context)\n",
    "\n",
    "plot(context=context)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
